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Abstract 

Folding of the triangular lattice in a discrete three-dimensional space is investigated by means of 
the transfer-matrix method. This model was introduced by Bowick and co-workers as a discretized 
version of the polymerized membrane in thermal equilibrium. The folding rule (constraint) is 
incompatible with the periodic-boundary condition, and the simulation has been made under the 
open-boundary condition. In this paper, we propose a modified constraint, which is compatible 
with the periodic-boundary condition; technically, the restoration of translational invariance leads 
to a substantial reduction of the transfer-matrix size. Treating the cluster sizes L < 7, we analyze 
the singularities of the crumpling transitions for a wide range of the bending rigidity K. We observe 
a series of the crumpling transitions at K = 0.206(2), —0.32(1), and —0.76(10). At each transition 
point, we estimate the latent heat as Q = 0.356(30), 0.08(3), and 0.05(5), respectively. 
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I. INTRODUCTION 



At sufficiently low temperatures, the polymerized membrane becomes flattened macro- 
scopically [jj; see Refs. Q-Q for a review. (The constituent molecules of the polymerized 
membrane have a fixed connectivity, and the in-plane strain is subjected to finite shear 
moduli. This character is contrastive to that of the fluid membrane [5|, la], which does not 
support a shear.) The flat phase is characterized by the long-range orientational order of 
the surface normals. It is rather peculiar that such a continuous (rotational) symmetry is 
broken spontaneously for such a two-dimensional manifold. To clarify this issue, a good deal 
of theoretical analyses have been reported so far. However, it is still unclear whether the 
transition is critical [7 19 1 or belongs to a discontinuous one with an appreciable latent heat 
-|22|. Actually, in numerical simulations, it is not quite obvious to rule out the possibility 



20 



of a weak-first-order transition 



23 



24] ; see also Ref. [25 ]. 

Meanwhile, a discretized version of the polymerized membrane was formulated by Bowick 
and coworkers 



26 



. To be specific, they considered a sheet of the triangular lattice em- 



bedded in a discretized three-dimensional space (face-centered-cubic lattice); see Fig. E] fa . 
(Even more simplified folding model, the so-called planar folding, was studied in Refs. 29- 
3JJ.) Owing to the discretization, the folding model admits an Ising-spin representation, for 



which a variety of techniques, such as the mean-field theory and the transfer-matrix method, 
are applicable. A peculiarity of this Ising magnet is that the spin variables are subjected 
to a local constraint (folding rule), which is incompatible with the periodic-boundary con- 
dition. Because of this difficulty, the open-boundary condition has been implemented so far 



26 



28 



32 



were considered 



34 



33l | . With the full diagonalization method, the finite clusters with the sizes L < 6 
26j, |28|. By means of the density- matrix renormalization group (DMRG) 



351 ]. the clusters with L < 29 and 26 were treated in Refs. [32j and 33j, respectively. 



(The results are overviewed afterward.) 

In this paper, we modify the local constraint, aiming to implement the periodic-boundary 
condition, and restore the translation invariance. Technically, the restoration of translation 
invariance admits a substantial reduction of the transfer-matrix size. Taking the advantage, 
we treat the sizes up to L = 7, and analyze the singularities of crumpling transitions in 
detail. 

The cluster-variation method (CVM) (based on the single-hexagon approximation) re- 
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vealed a rich character of the discrete folding 27|, |28|. According to CVM, there appear the 
totally flat, octahedral, tetrahedral, and piled-up phases, as the bending rigidity K changes 
from oo to — oo. Namely, in respective phases, the triangular-lattice sheet crumples up to 
form a octahedron, a tetrahedron, and a triangular plaquette; see Fig. 3 of Ref . 26] . (This 
picture is based on a single-hexagon approximation of CVM. Beyond a mean-field level, the 
thermal undulations may be induced, particularly, in the vicinity of the transition point, 
disturbing the shape of the crumpled sheet significantly.) More specifically, the crumpling 
transitions separating these phases are estimated as K = 0.185, —0.294, and —0.852; here- 
after, we abbreviate the set of parameters as (0.185,-0.294,-0.852). At each transition 
point, the latent heat is estimated as (0.229, 0.14, 0); namely, the third transition is continu- 



ous according to CVM. On the one hand, the DMRG simulation [32|, [33J indicates the transi- 
tion point [0.195(2), -0.32(1), -0.76(1)] with the latent heat [0.365(5), 0.04(2), 0.03(2)]; the 
character of the third transition point is still controversial. 

The rest of this paper is organized as follows. In Sec. [Til we propose a modified folding 
rule (OQ); the transfer-matrix formalism {26 1 is explicated in the Appendix. In Sec. IHIl we 
present the numerical results. The singularities of the crumpling transitions are analyzed in 
detail. In Sec. IIV) we present the summary and discussions. 



II. 



A MODIFICATION OF THE FOLDING RULE 



In this section, we present a modified folding rule, Eq. (JT]). As mentioned in the Introduc- 
tion, the folding rule, enforced by the prefactors in Eq. flA.lj) . is too restrictive to adopt the 
periodic-boundary condition. So far, the numerical simulation has been performed under 



the open-boundary condition 



26 



28 



32 



33). 



To begin with, we outline the transfer-matrix formalism; an explicit algorithm is presented 
in the Appendix. According to Ref. jiffl ]. through a dual transformation, the triangular- 
lattice folding reduces to an Ising model on the hexagonal lattice, Fig. [9] (a). A drawing of 
a transfer-matrix strip is presented in Fig. [9] (b). A peculiarity of this reduced Ising model 
is that the spins surrounding each hexagon are subjected to a constraint. (The constraint 
originates from the folding rule.) To be specific, the prefactors UjVj(= 0, 1) of the transfer- 
matrix element flA.ip restrict the configuration space. As mentioned in the Introduction, this 
constraint is incompatible with the periodic-boundary condition. (For instance, a cylindrical 
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paper supports a large strain, whereas an open paper is flexible.) 

Aiming to restore the translational invariance, we modify the prefactors. We replace Eq. 
(jATT]) with 

1 L 

1 Etfl U i V i) [(1 - P)UM + p] exp[- Hk{K) - H^K')]. (1) 

i=i j^i k^i 

Here, the parameter L denotes the system size, and the explicit formulas for the constraint 
UjVj and the elastic energy Hk(K) are shown in the Appendix. As compared with the 
original form flA.ll) . our modified expression ([1]) has a defect at I, where the folding-rule 
constraint is released. Because the defect is distributing uniformly by the summation Y^i=i 
and the normalization factor 1/L, the translational invariance is maintained. The param- 
eter K' describes the local elastic constant at the defect. The probability of the defect is 
controlled by the parameter p; at p = 0, the original constraint V\V\ recovers, whereas at 
p — 1, the constraint disappears XJ\V\ — > 1. We stress that a single defect does not alter the 
thermodynamic (bulk) properties. As a byproduct, two tunable parameters p and K' are 
available. The parameters are adjusted to 

(p,K') = (0.7, 1.5K). (2) 

A justification of this choice is given in Sec. IIII Bl (Fig. 

III. NUMERICAL RESULTS 

In this section, we present the simulation results. We employed the transfer-matrix 
method (Appendix) with a modified folding rule, Eq. ([T]). The numerical diagonaliza- 
tion was performed within a subspace specified by the wave number k = and the parity 
even. (In a preliminary survey, we confirmed that the dominant-eigenvalue (thermal equi- 
librium) state belongs to this subspace.) Here, we make use of the spin- inversion symmetry 
(Ti,Zi — > ±cTj,±2j. (This symmetry group originates from the overall rotation of the crum- 
pled triangular-lattice sheet.) We stress that the wave number k makes sense owing to the 
restoration of the translational invariance (JT]). 
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A. Crumpling transitions: A preliminary survey 



In Fig. [U we plot the free-energy gap 

A/ = h ~ fi, (3) 

for the bending rigidity K. Here, the free energy per unit cell is given by fi = — lnAj/(2L) 
with the (sub) dominant eigenvalue A^) of the transfer matrix. [Here, the unit cell stands 
for a triangle of the original lattice rather than a hexagon of the dual lattice; see Fig. (a).] 
As mentioned in the Introduction, the triangular-lattice sheet becomes crumpled, as the 
rigidity K changes from oo to — oo. In Fig. [TJ we see a number of signatures of the 
crumpling transitions around K pa 0.2, —0.3, and —0.8. (Note that the closure of the 
excitation gap indicates an onset of phase transition.) On the one hand, the CVM analysis 



271 . |28| predicts a series of crumpling transitions at (0.185,-0.294,-0.852). The results 
appear to be consistent with those of Fig. [JJ suggesting that the excitation-gap closure 
indicates a location of the crumpling transition. Detailed analyses of each singularity are 
made in Sec. IIII CI and Sec. IIIIDl 



B. Simulation at p = 1 and K' = 

As a comparison, we provide a simulation result, setting the defect parameters to p — 1 
and K' = tentatively. This parameter set has an interpretation that a rupture (pair of open 
edges) distributes uniformly along the transfer-matrix strip. (This situation is an extention 
of the open-boundary condition, for which the rupture is static.) 

In Fig. El we present the free-energy gap Af for the bending rigidity K; the range of K 
is the same as that of Fig. [TJ The signatures of crumpling transitions in Fig. [2] are less clear, 
as compared to those of Fig. [TJ This result indicates that the choice of the defect parameters 
affects the finite-size behavior. In the preliminary stage, we surveyed a parameter space of p 
and K', and arrived at a conclusion that the above choice, Eq. (j2j), is an optimal one. Note 
that these parameters are the byproduct of the modification of the folding rule, Eq. ([TJ). 
Here, we make use of these redundant parameters so as to improve the finite-size behavior, 
aiming to take the thermodynamic limit reliably. 
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C. Crumpling transition in K > 0: Analysis of the latent heat via the Hamer 
method 



In Sec. IIII A[ we observed a series of crumpling transitions. In this section, we analyze 
the singularity of a transition in the K > side. 

To begin with, we determine the transition point. In Fig. |3l an approximate transition 
point K C (L) is plotted for 1/L 2 . Here, the approximate transition point minimizes the 
free-energy gap 

d K Af\ K=Kc{L) = 0. (4) 

The least-squares fit yields an estimate K c = 0.20617(70). Similarly, as for 5 < L < 7, 
we obtain an estimate K c = 0.20506(99). A discrepancy between these results may be an 
indicator of possible systematic error. The systematic error appears to be comparable to 
the fitting error. Regarding both errors as the sources of error margin, we obtain 

K c = 0.206(2). (5) 

Based on the transition point K C (L), we estimate the amount of latent heat. According 
to Ref. |36| |. the low-lying eigenvectors of the transfer matrix contain information on the 
latent heat. We explain the underlying idea, and present the scheme explicitly. At the 
discontinuous (first-order) transition point, the low- lying spectrum of the transfer matrix 
exhibits a level crossing, and the discontinuity (sudden drop) of the slope reflects a release 
of the latent heat. However, the finite-size artifact (level repulsion) smears out the singu- 



larity. According to Ref. [36j], regarding the low-lying levels as nearly degenerate, one can 
resort to the perturbation theory of the degenerated case, and calculate the level splitting 
(discontinuity of slope) explicitly. To be specific, we consider the matrix 

vJ V » V A, (6) 

\ v 2 i v 22 J 

with Vij = (i\dj(T\j). The matrix T denotes the transfer matrix; namely, the matrix element 
dxT is given by a product of Eq. ( lA.lj) and Eq. ( 1A.5I) with K dropped. The bases, |1) 



and |2), are the (nearly degenerate) eigenvectors of T with the eigenvalues A 12 , respectively. 
The states {\i)} are normalized so as to satisfy (i\T\i) = 1. According to the perturbation 
theory, the eigenvalues of Eq. fl6]) yield the level-splitting slopes due to K. Hence, the latent 
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heat (per unit cell) is given by a product of this discontinuity and the elastic constant 

Q(L) = \K c (L)\y/(V n ~ V 22 ) 2 + AV 12 V 21 ^ (7) 

for the system size L. 

In Fig. HI we plot the latent heat ([7]) for 1/L 2 . The least-squares (linear) fit yields an 
estimate Q = 0.356(12) in the thermodynamic limit. Similarly, we obtain Q = 0.3774(93) 
for 5 < L < 7. Again, the systematic error appears to be comparable to the fitting error. 
We estimate 

Q = 0.356(30). (8) 
This is a good position to address a few remarks. First, the latent heat, Eq. ([8]), agrees 



with that of DMRG 



321 ] (see the Introduction), whereas the transition point, Eq. 



lies out of the error margin. This discrepancy may indicate an existence of systematic 
error as to the determination of K c . A peculiarity |30[ of this transition is that in the 
K > K c side, the system becomes completely flattened; namely, there exist no thermal 
undulations, as if the system is in the low-temperature limit K — > oo. This peculiarity may 
gives rise to a bias to K c . (As a matter of fact, in Ref. [32| . a pronounced hysteresis was 
observed.) On the contrary, we confirmed that the ambiguity of K c does not influence the 
latent heat Q very much. (Because of this K c independence, the result Q is reliable, as the 
above-mentioned consistency with DMRG suggests.) In the next section, the remaining two 
transitions are analyzed in a unified manner. Second, we consider the 1/L 2 extrapolation 
scheme. The finite-size data converge rapidly to the thermodynamic limit around the first- 
order transition point, because the correlation length (typical length scale) £ remains finite. 
Hence, the dominant system-size corrections should be described by 1/L 2 (rather than 1/L). 



D. Crumpling transitions in K < 

In this section, we analyze the remaining transitions in the K < side. 

First, we analyze the transition around K pa —0.3. In Fig. [5j we plot the transition 
point (jlj) for 1/L 2 . The least-squares fit to these data yields an estimate K c = —0.320(12). 
Similarly, we obtain K c = —0.316(28) for 5 < L < 7. The systematic error appears to be 
negligible, as compared to the fitting error. Considering the latter as the source of error 
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margin, we estimate the transition point as 



K c = -0.32(1). (9) 

In Fig. [6j we plot the latent heat © for 1/L 2 . The least-squares fit yields Q = 0.077(15). 
Similarly, for 5 < L < 7, we obtain Q = 0.058(28). The fitting and systematic errors are 
comparable. Considering them as the sources of error margin, we estimate the latent heat 
as 

Q = 0.08(3). (10) 

Second, we turn to the analysis of the transition around K w —0.8. In Fig. [71 we plot 
the transition point (J4j) for 1/L 2 . The least-squares fit to these data yields an estimate 
K c = -0.76(10). Similarly, we obtain K c = -0.72(24) for 5 < L < 7. The fitting error 
dominates the systematic error. Neglecting the latter, we obtain 

K c = -0.76(10). (11) 

In Fig. El we plot the latent heat © for 1/L 2 . The least-squares fit yields Q = 0.049(51). 
Similarly, for 5 < L < 7, we obtain Q = 0.01(11). Again, the systematic error appears to 
be negligible. We estimate the latent heat as 

Q = 0.05(5). (12) 

Last, we consider a shaky character of Figs. [3HH] (in particular, Figs. [MS)- Such a shaky 
character is an artifact due to the cluster size. In the preliminary stage, we surveyed the 



planar folding 29h31| for considerably large system sizes L < 14; the configuration space 
of the planar folding is much restricted. As a result, we found that the finite-size behavior 
is irregular with respect to L; this irregularity is an obstacle to making an extrapolation. 
Roughly speaking, the finite-size behavior is categorized by L = 0, 1, and 2 mod 3. Although 
the enlarged configuration space suppresses this irregularity, a slight irregularity for L = 6(= 
mod 3) still remains. A slight bump of Figs. [5]and[7]may be due to this irregularity. Hence, 
we consider that the deviations (seemingly curved plots) in Figs. |3]|8]are not systematic ones. 
Rather, considering them as a source of errors, we estimate the error margin by making two 
independent extrapolations for different sets of system sizes. 
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IV. SUMMARY AND DISCUSSIONS 



We proposed a modified folding rule, Eq. flTJ, which enables us to simulate the discrete 
folding without the open edges. By means of the transfer-matrix method, we investigated 
a series of the crumpling transitions. We estimate the transition point and the latent heat 
as [0.206(2), -0.32(1), -0.76(10)] and [0.356(30), 0.08(3), 0.05(5)], respectively. Our result 
agrees with the preceding CVM and DMRG results. In particular, our result is in quantita- 
tive agreement with that of DMRG. 

According to Ref. j^, the third singularity (around K m —0.8) is so subtle that it 
could not be captured until L = 8. On the contrary, our data in Figs. [TJ and [7] succeeds 
in detecting a signature of a crumpling transition even for 4 < L < 7. As anticipated, the 
restoration of the translation invariance leads to an improvement of the finite-size behavior. 

As mentioned in the Introduction, it is still unclear whether the third transition is con- 
tinuous or belongs to a weak- first-order transition [33] . The present result (Tl2|) does not 
exclude a possibility of a continuous transition. According to Ref. [28j, around K « —0.8, 



through a truncation of the configuration space, the discrete 
version of the folding model, the so-called planar folding 29 



biding reduces to a simplified 



31], which exhibits a continu- 



ous transition in the K < side. An examination of this truncation process may provide 
valuable information on the nature of this phase transition. This problem will be addressed 
in the future study. 



Appendix: Transfer-matrix formalism for the discrete folding 



In this Appendix, we present the transfer-matrix formalism for the discrete folding [26]. 
Before commencing a mathematical description, we explicate a basic feature of the discrete 
folding. We consider a sheet of the triangular lattice, Fig. [9] (a). Along the edges, the sheet 
folds up. The fold angle 9 along the edges is discretized into four possibilities, namely, "no 
fold" (9 = tt), "complete fold" (9 = 0), "acute fold" [9 = arccos(l/3)], and "obtuse fold" 
[9 = arccos(— 1/3)]; in other words, the triangular lattice is embedded in the face-centered- 
cubic lattice [261 ] . 

The above discretization leads to an Ising-spin representation. The mapping, the so- 



called gauge rule, reads as follows [26|. We place two types of Ising variables {t7j,Zj} at 
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each triangle i (rather than each joint); see Fig. [9] (a). Hence, hereafter, we consider the 
a spin model on the dual (hexagonal) lattice. The gauge rule sets the joint angle between 
the adjacent triangles. That is, provided that the z spins are antiparallel [z\Z 2 = -1) for 
a pair of adjacent neighbors, the joint angle is either an acute or obtuse fold. Similarly, 
if 0\0 2 = —1 folds, the relative angle is either a complete or obtuse fold. The spins are 
subjected to a constraint (folding rule); The prefactors UjVj of the transfer-matrix element 
flA.lj) enforces the constraint. 

As a consequence, the discrete folding reduces to a two-component Ising model on the 
hexagon lattice. Hence, the transfer-matrix strip looks like that drawn in Fig. M (b). The 
row-to-row statistical weight T^ a . tZ .y^ cr / tZ q yields the transfer-matrix element. The transfer- 
matrix element for the strip width L is given by 26 ] 



with 



and 



L 

T {z ^ Uz> > } = ([[ UM) eM-H/T), (A.1) 



Uj = 5(cr 2 j-2 + &2j-i + cr 2j + (y' 2j -i + cr' 2j + cr' 2j+1 mod 3, 0), (A. 2) 



Vj = Yl S(a c (z 2 j, zjy-i, 4i-u Ay z 2j+i) mod 2, 0). (A.3) 

c=l 

The factors {Uj,Vj} enforce the constraint (folding rule) as to the spins surrounding each 
hexagon. Here, 5(m,n) denotes Kronecker's symbol, and a c is given by 

6 ^ 

a c (^i, . . . ,z 6 ) = ^-(1 - ZiZ i+1 )5C^2aj mod 3,0). (A.4) 

i=l j=l 

The Boltzmann factor exp(—H/T) is due to the bending-energy cost H. Hereafter, we 
choose the temperature T as a unit of energy; namely, we set T = 1. As usual, the bending 
energy is given by the inner product cos^- of the surface normals of adjacent triangles. 
Hence, the bending energy is given by a compact formula 

H = Y J H k {K)- (A.5) 

k 

Here, the index k specifies each hexagon as in Eq. (lA.ip . The local energy of each hexagon 
iffc is given by 

6 6 j 

H k (K) = -0.5 J2 K cos ki+i = -0.5 J]-X(7 i (7 i+1 (l + 2 %+ i), (A.6) 

i=l i=l 
10 



with the bending rigidity K. Here, the summation X^=i runs over an ver tices around the 
hexagon k. (The overall factor 0.5 compensates the duplicated sum.) 
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FIG. 2: The free-energy gap ([3]) is plotted for the bending rigidity K and the system size L =(+) 
4, (x) 5, (*) 6, and (□) 7. Tentatively, the defect parameters [Eq. ([1])] are set to p = 1 and K' = 0. 
The data appear to be scattered, as compared to those of Fig. [TJ This result indicates that the 
finite-size behavior is improved by adjusting the defect parameters. 
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FIG. 3: The crumpling transition K ~ 0.2 observed in Fig. [I] is analyzed in detail. The transition 
point K C (L) © is plotted for 1/L 2 . The least -squares fit yields an estimate K c = 0.20617(70). 
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FIG. 4: The latent heat Q(L) © is plotted for 1/L 2 . The least -squares fit yields an estimate 
Q = 0.356(12). 
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FIG. 5: The crumpling transition K ~ —0.3 observed in Fig. [T]is analyzed in detail. The transition 
point K C (L) flU) is plotted for 1/L 2 . The least -squares fit yields an estimate K c = —0.320(12). 
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FIG. 6: The latent heat Q(L) © is plotted for 1/L 2 . The least -squares fit yields an estimate 
Q = 0.077(15). 
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FIG. 7: The crumpling transition K ~ —0.8 observed in Fig. [T]is analyzed in detail. The transition 
point K C (L) (jU) is plotted for 1/L 2 . The least-squares fit yields an estimate K c = —0.76(10). 
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FIG. 8: The latent heat Q(L) © is plotted for 1/L 2 . The least -squares fit yields an estimate 
Q = 0.049(51). 
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(a) 




(b) 



FIG. 9: (a) We consider a discrete folding of the triangular lattice. In order to specify the fold 
angle, we place two types of Ising variables {#j,Oi} on each triangle i rather than at each joint 



(gauge rule 26(] ) . Hence, hereafter, we consider a spin model on the dual (hexagonal) lattice, (b) 
Construction of the transfer matrix. The row-to-row statistical weight yields the transfer-matrix 
element, Eq. (|A.ip . So far, the open-boundary condition has been imposed. Here, we restore the 
translational invariance by using a modified folding rule ([T]). 
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